February 2, 2008 



A modified Balitsky-Kovchegov equation 



E. Gotsman, * fl ) E. Levin,t U. Maoi 4 a ) and E. Naftali § a ) 



a ) HEP Department 
School of Physics and Astronomy 
Raymond and Beverly Sackler Faculty of Exact Science 
Tel Aviv University, Tel Aviv, 69978, ISRAEL 

5 ) DESY Theory Group 
22603, Hamburg, GERMANY 



February 2, 2008 



Abstract 

We propose a modified version of the Balitsky-Kovchegov (B-K) evolution equation, 
which includes the main NLO corrections. We use the result that the main NLO corrections 
to the BFKL kernel are the LO DGLAP corrections. We present a numerical solution of the 
modified nonlinear equation, and compare with the solution of the unmodified B-K equation. 
We show that the saturation momentum has a sharp increase in the LHC energy range. Our 
numerical solution shows that the influence of the pre-asymptotic corrections, related to 
the full anomous dimensions of the DGLAP equation, are rather large. These corrections 
moderate the energy behaviour of the amplitude, as well as the value of the saturation scale. 
All our calculations are made with a fixed value of a,. 
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1 Introduction 



High density QCD has a long history, starting with Refs. P3I2HH], and is now entering a new phase 
of its development: direct comparison with the experimental data This makes stringent 

demands on the theoretical approach to provide reliable predictions. The main theoretical tool is 
the non-linear Balitsky-Kovchegov (BK) evolution equation |Hj which has a few deficiencies: 

• It is correct only in the leading log(l/x) approximation of perturbative QCD (pQCD). In 
practical terms it means that the kernel of the BK equation should be the BFKL kernel 
in the leading order (LO) of pQCD; 

• This equation is the mean field approximation to the JIMWLK equation jH], which is more 
general, but also more complicated. The JIMWLK equation has not been sufficiently well 
investigated theoretically, to attempt any description of the experimental data; 

• We know that the BK equation is not correct in the saturation regional HUj. hence cannot 
be viewed reliable tool to explore this region; 

• The region where we can neglect non-linear corrections, should be specified by conditions 
which are beyond the BK equation jTTj. 

In this paper we address only one of the above problems, namely, the higher order corrections to 
the kernel of the BK equation. 

As it has been discussed [U 12 El El , the non-linear corrections should be incorporated before 
taking into account the next-to-leading order (NLO) corrections to the kernel of the linear equation 
(BFKL or DGLAP ^2]). The non-linear corrections are proportional to the exchange of two BFKL 
ladders [H E] . Since the dipole amplitude in the BFKL approach has the following behaviour 

M nr n/ 4 Q 2 ^BFKL fl T\ 

1 ' non-linear term ' \ J 

where, N denotes the dipole-target amplitude and behaves as 

-^linear term OC Og S^ BFKL (1-2) 

and Abfkl is the intercept of the BFKL Pomeron (Abfkl oc as in the leading order). 

Comparing Eq. (jl.ljl and Eq. (jl.2j) one sees that for a wide range of energies one has to include 
the non-linear effects when 

y = Ins > — In ( — ) (1.3) 
a s \a s J 

On the other hand, the corrections of the order of a s in Abfkl, becomes essential only for 
energies a s y ~ 1 or ?/ = Ins > -4-. Comparing this estimate with Eq. (jl.3|) we see that the 
correct theoretical strategy is to take into account all non-linear corrections first, and then to 
calculate the NLO corrections to the BFKL kernel. 
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However, explicit calculations of the NLO corrections to the BFKL kernel ^3], show that a| 
corrections to A BFKL are rather large, and for the past five years, have been the subject of 
many debates and discussions ^3]. The correct resummation of the NLO corrections have been 
performed by the Florence group (see Ref. j201), and it turns out that the NLO corrections are 
essential in estimates of the non-linear effects. As has been discussed in Refs. [TJ UHl EH HZ| the 
estimates for the saturation momentum does not depend on the precise structure of the non-linear 
corrections. Mueller and Triantafyllopoulos found in Ref. [HI] > that the NLO corrections crucially 
change the energy behaviour of the saturation scale, making it close to the phenomenological 
saturation scale appearing in the Golec-Biernat and Wusthoff model 

At first sight, it seems that we are in a vicious circle: on the one hand, the NLO corrections 
are important, while on the other hand, they should be taken into account only after non- linear 
corrections. At the moment, we have no idea of how to combine these two effects into one 
equation. In this situation both the non-linear BK equation in the LO, and the linear BFKL 
equation in NLO are unreliable. 

In this paper we suggest a way out of this dilemma. We propose a modified version of the BK 
equation that takes into account the main NLO corrections. The idea on which this equation is 
based, is the fact that the main next-to-leading order corrections to the BFKL kernel are actually 
the leading order DGLAP corrections [201 Exj ■ Therefore, in our equation we suggest a procedure 
in which the non-linear evolution equation based on the LO BFKL kernel, is combined with LO 
DGLAP evolution. Different approaches of how to make such matching, have been discussed 
both for the linear equation (see for example Refs. PIE] ), and for the non-linear equation (see 
Refs. |22*1 123} l2l| 125] and references therein). We trust that we have found an economical and 
simple method to include LO DGLAP corrections in the framework of the non-linear equation. 
The equation is discussed in detail in section 2. Section 3 is devoted to the numerical solution 
of the modified non-linear equation, and to a comparison with the solution of the BK equation 
without any modification. In the Conclusion section 4, we summarize the results and discuss the 
application of the new equation to a global fit of the experimental data. 



2 The modified non-linear equation 
2.1 The saturation scale 

In this section we discuss the value of the saturation scale, which does not depend on the form of 
the non-linear term in the BK equation. The linear part of the evolution equation can be written 
in a simple form, if we go to the double Mellin transform for the dipole scattering amplitude (N), 
namely, 

N(t = Hr 2 /R 2 ),Y = \n(l/x) ] b) = [ ±L N ^ r b) e »r + 7« (2 .l) 

where r and R are, respectively, the sizes of the dipole and the target, x denotes the Bjorken 
variable for the dipole-target scattering, and b is the impact parameter for the reaction. The 
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BFKL equation (or in other words the linear part of the BK equation) has the form: 

u = u(y) = asXLoil) + a 2 s XNLo(l)- (2.2) 
For Xlo(i) we have the well known expression jjj: 

Xlo(i) = 2^(1) - ^(7) - V(l"7) = ^ + 73^ + X ^(7), (2-3) 

with ip = dlnr(7)/d7 and T the Euler Gamma function. We will discuss the form of Xnlo(i) 
later. Xlo{i) = 2^(1) — if)(l + 7) — if) (2 — 7) denotes the contribution of the higher twist 
which we will discuss below. 

The form of Eq. ()2.2|) is clear since the l.h.s. and r.h.s. of the BFKL equation 

dN(Y,0 



BY 

have the following Mellin transforms: 

dN(Y,£) 
dY 



j K BFKL (Z,ON(Y,Odt (2.4) 



uN(w,i) (2.5) 



f K BFKL ^ON(Y,C)di - u(n)N{u>,i). (2.6) 

The saturation momentum can be calculated without making any assumption regarding the char- 
acter of the non-linear corrections, and it has the following form: 



Q 2 S (Y) = Q 2 s (Y )e W (^-(Y-Y ) - , w| ' , hMY/Y, 

\ J- 7cr 



2tt , 1 



(1 -7cr) 2 V W"(7cr) WY 




(2.7) 



where Y = ln(l/x) is our energy variable, ^"(7) = d 2 uj('-f)/(d'j) 2 , the value of 7 cr can be found 
from the equation [H IT3]: 

1 7cr *^7cr 

we have normalized the value of the saturation momentum at Y = Y . 

The first term was given in the GLR paper pp, the second in Ref. and the third in Ref. 

[15]. Fixing the value of Q S (Y ) = 0.37 GeV 2 at x = 10 _2 (F = 4.6), we obtain the saturation 
momentum shown in Fig. ^for a$ = 0.2 (see curve 1). We see that the saturation momentum 
shows a steep rise towards the LHC energies. Consequently, prior to any discussion of the value 
and energy dependence of the dipole scattering amplitude, we need to understand how the next- 
to-leading order corrections could effect the value of the saturation scale. As we have mentioned, 
these corrections change the value of the scale considerably [T§1 IT7] . This can be seen in Fig. Q 
comparing curve 1 with curve 3 which shows the behaviour of the saturation scale for the NLO 
BFKL kernel. 
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2.2 The NLO BFKL kernel 



We now discuss the NLO corrections, and suggest a procedure to include these in the equation, 
in a manner so that collective effects are treated properly. The NLO BFKL kernel is known [13j, 
and all theoretical problems related to this kernel have been solved [20J. 

Our procedure is based upon the important observation, made by the Durham group [Hj, ac- 
cording to which the main NLO corrections to the BFKL kernel, can be taken into account using 
a simple expression for the NLO BFKL kernel: 

, s fl+ujAAuj) 1 1 + ujAAuj) 1 \ HT/ . . . 

asXNLoil) = ^-- + - jr 1 --, -wXioT, 2.9 

\ 7 7 1 — 7 + c<j 1 — 7 j 

where, 

A(u) = -11/12 + 0(u) + n F (^§- ?, G ( W ) P gG (o;) - 1/3^ , (2.10) 
with P(u;) being the DGLAP kernel. 

The singularities in Eq. (|2.9|) describe the different branches of evolution corresponding to the sizes 
of the interacting dipole (or the transverse momenta of partons). The pole at 7 = corresponds 
to the normal twist-2 DGLAP contribution, with the ordering in the transverse parton momenta 
Q > ...ki >t > k t)i+ i > ••• > Q , where Q is the typical virtuality of the target Q ~ 1/R. 
The pole at 7 = 1 in Eq. (|2.9|) corresponds to inverse k t ordering {k i)t < k tji+ i < . . . Q ). The 
other poles, at 7 = —1, —2, ... (7 = 2, 3, ... ), are the higher twists contributions due to the gluon 
reggeization 1 . 

As can be seen in Eq. (J2.9)) . the main alterations that we need to make to incorporate NLO, 
are in the sector which corresponds to the ordinary DGLAP evolution. We need to change the 
BFKL anomalous dimension, so as to account for the DGLAP anomalous dimension. The second 
important change that the NLO calculations induce, is in the inverse evolution. However, this 
branch of evolution is moderated by the non-linear effect, rather than by the NLO corrections. 
Indeed, it was shown in Ref. j2E] that the term 1/(1—7) leads to exponentially small corrections 
in the saturation region, or, in other words, for 7 > 7 cr « 0.37. Therefore, in attempting to find 
a new kernel for the non-linear equation, we can neglect changes in the inverse evolution, and 
keep the BFKL kernel without a shift from 7 to 7 — uj. 

Using the observation of Ref. [2Z], according to which Eq. (J2.10|) for Ai(u) can be well approxi- 
mated to an accuracy (> 95%), by Ai(u) = 1. We rewrite the NLO BFKL kernel in a very simple 
form: 

Xnlo(i) = -uxlo{i)- (2.11) 

It should be mentioned that A\{u) = 1 corresponds to the expression for the anomalous dimension 
of the DGLAP equation: 



7loM = 




(2.12) 



1 Onc of us (E.L.) acknowledges long and fruitful discussions with J. Bartels on the physical meaning of all 
parts of the BFKL kernel. 
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which approximates the DGLAP anomalous dimension in leading order to within an accuracy of 
95%. The actual deviation of Eq. (J2.12)) from the correct expression for the DGLAP anomalous 
dimension in leading order, is even less than 5%. 

The saturation momentum for Eq. (|2.11|) is plotted in Fig. ^ (see curve 2). This momentum 
is larger than the saturation momentum given by Eq. (|2.9j) . but we trust that the non- linear 
corrections will suppress it. 

2.3 Modified non-linear equation 




i i i i i i i i i i i i i i i i i i i i iJ 

6 8 10 12 14 

Y 

Figure 1: The saturation momentum as a function of x (Y = ln(l/x)j in the LO BFKL (curve 1) , for 
Eq. \2.11\) (curve 2) and for Eq. \2. ,9J) (curve 3) for as = 0.2. Curve 4 is the Golec-Biernat 
and Wusthoff phenomenological saturation momentum. 

Using Eq. 1)2.11)1 we obtain the full kernel for the linear equation in the form: 

w(7) = ®sXlo(i) (1 - w) (2.13) 
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This kernel imposes energy conservation (see Ref. |2Z])> and describes the NLO BFKL kernel. It 
does not include the contribution coming from inverse ordering, which should be suppressed in 
the solution to the non-linear equation. 

Our suggestion is to use Eq. (|2.13J) as the kernel for the non-linear equation. 
The modified non-linear equation with this kernel has the form 

ON (r, Y; b) C F a s f d 2 r> r 2 ON (r',Y;b - \(f- f')^ Cpas r d 2 r > r 2 



(^2N (^r',Y;b - ~ (r - r ')J - JV (V, K; 6 - ~ (r - f ') J N r K; 6 - ~f 'j j . (2.14) 

Using Eq. (|2.5j) and Eq. (|2.6j) one can recognize , that the r.h.s. of the equation is the usual 
BK equation, while the second term on the l.h.s. of the equation is a new contribution which 
takes into account our modification of the BFKL kernel. Indeed, the term a$ & Xlo{i) h as the 
following form in the Y, r representation 

a s u X Lo{l) - a s J K LO (Y,r')d 2 r'^^l. (2.15) 

In Eq. (I2.14j) the linear part conserves energy, however, in the non-linear part energy is not con- 
served. It should be stressed that although Eq. ()2.14j) is written in the leading order of as, it also 
includes terms that are not proportional to (asln(l/x)) n . Such resummation is legitimate, and 
this type of correction could be taken into account before considering the non-linear corrections. 

In the dipole approach the kernel of the non-linear equation is the probability for one dipole 
(with size r) to decay into two dipoles (with sizes r' and \r — r'\). Assuming that this physical 
interpretation is correct in the NLO, we can rewrite Eq. (J2.14|) in a more general form: 



dN{r,Y;b) C F a s f d 2 r' r 2 



Qy 7T 2 J {f- r') 2 r' 2 

2Nlr',Y;b- ^(r-f'^j - N (r',Y; b - i (r - r ') j JV (f - f ',Y; b - ~f ' 

° 2N ( r' Y:b - - (r- f')] - JV \r\Y\ b— - (f- r') \ JV If- f',Y: b — -r ' 



(2.16) 



dY V V 2 v 7 V 2 v 7 V 2 

The advantage of this equation is that energy is conserved by the non-linear term, as well as by 
the linear one. 
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2.4 Pre-asymptotic corrections to the value of the saturation scale 



Before solving Eq. (|2.1fi|) numerically, we would like to justify the need for a numerical solution. 
Indeed, at first sight we have Eq. ()2.7|) (see also Eq. ()2.17|) below) for the saturation scale, and 
we have a good understanding of (i) how a solution to BK equation approaches the saturation 
boundary N — 1 (see for example Ref. [2^] ), and (ii) how the amplitude behaves in the vicinity 
of the saturation scale [TBI EH] ■ 

We believe that the main shortcoming of our analytical approach at present, is the fact that it can 
only be applied for asymptotically high energies, while for any practical use we need to know the 
solution at rather low energies. To specify what high and low are (for this particular problem), 
we return to the discussion of the value of the saturation scale given by Eq. ()2.7|) . However, we 
now consider the equation from a different point-of-view, and ask ourselves, how large are the low 
energy corrections to the value of the saturation scale. 

Fig. |21 shows the saturation scale when only the high energy contribution (the first term in 
Eq. ()2.7|) ) is taken into account (curve 1), and the saturation scale including the low and high 
energy corrections (curve 2). One can see that the difference is very large, even for the LHC 
energy range. The same situation holds for the case of the running as- In this case the formula 
for the saturation scale with the low energy corrections is: 

Ql(Y) = (2.17) 

6 ^ 4 \^fblcrX(lcr) J J 

The energy behaviour of the first term of Eq. (|2.17|) (high energy), and for both terms are shown 
in curves 3 and 4 of Fig. |2J The difference between these two curves, which is large, indicates the 
size of the low energy correction. 

The same conclusion regarding the essential low energy corrections to the value of the saturation, 
and its energy behaviour, has been derived from the numerical solution of the BK equation. In 
Ref. the energy dependence of the saturation scale in the LHC range of energies turns out 
to be more moderate, than predicted by asymptotic formulae of the type of Eq. ()2.7j) . 

We can therefore conclude, that to provide reliable predictions for the LHC range of energies, it 
is necessary to solve the BK modified equation numerically at low energies. 

The second conclusion, which we can derive from Fig. comparing curve 4 and curve 2, is 
the fact that the inclusion of the running QCD coupling is essential, so as to provide a reliable 
extrapolation of the HERA data to the LHC energy range. 

3 Numerical Solution 

Ideally, our primary goal would have been to obtain a numerical solution to the full, energy 
conserving, modified BK equation. Although, Eq. (j2.16j) is a well defined integro-differential 
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Figure 2: The saturation momentum as function of x (Y = ln(l/x) ). Saturation momenta for the 
case of fixed QCD coupling (as = 0.2): curve 1 shows the high energy behaviour of Q s (the 
first term in Eq. \2. 7| ) ), while curve 2 shows all term in Eq. \2.T\ l. Curves 3 and 4 are 
the same as curve 1 and 2 but for running as- Curve 3 shows the first term of Eq. \2.1T\j 
which describes the high energy behaviour of the saturation scale curve 4 is calculated using 
all terms of Eq. {Offi . 

mathematical equation, obtaining a numerical solution to this equation is a rather complicated 
task, due to the appearance of ON/dY on both sides of the equation. 

A natural first step towards a numerical solution of Eq. (|2.1fijl is to rewrite the DGLAP correction 
in the double log approximation (the second term in l.h.s. of Eq. ()2.14)l ). including only leading 
twists, and neglecting the impact parameter dependence. We follow Ref. [2E]> an d modify the 
kernel of the BK terms in the modified equation so as to exclude contributions of inverse fan 
diagrams. 
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We write the evolution equation in a differential form as: 



dN(r,Y;b) C F a s 



dY 



2r z 



7T- 



d 2 r' ON (r', Y; b) 



„/4 



dY 



+ 



d 2 r' r 2 



B(R-r')Q(R- \r-r'\) x 



2N (r',Y;b-\{r-?)\ - N (r',Y;b- ±(r-P)) N (r - ?,Y;b - \A \. (3.1) 



In our numerical calculations we have taken the value of the target size, R, to \/l0 GeV -1 . Note 
that, in principle, in the double log approximation, DGLAP contributions should be calculated 
for large distances, irrespective of the size of the target. However, as our pure BK solution can 
only be trusted for r < R (see j2Hj for further details), we impose the same restriction on the r' 
integration of the first term of the r.h.s of Eq. (j3.1|) . 

Our solution to Eq. (|3.1|) will be presented in terms the integrated quantity: 



N(Y,r;b) = df\\ N(Y,r;b;fu) 



(3.2) 



where r = \r\, b = |6|, and rn = b ■ r/(br). This quantity determines the physical observables in 
DIS, such as F 2 and gluon densities. 

Briefly, our numerical procedure is as follows. Denoting a particular rapidity by , at which the 
solution and its derivative are known, we define the solution at Y i+ i = Y^ + h(Y) as a matrix, in 
which the matrix elements correspond to (fixed Y) solutions at different r and b, integrated over 
f\\: 

dN(Yr,',b',r\\) 



N(y i+1 , r; b) = N(Xi, r; b) + h(Xi) / df\ 



dY 



(3.3) 



Y=Y l 



where dN(Y, r; b; f\\)/dY is given by the r.h.s. of ()3.1|) . and h(yi) is a variable step size in the 
rapidity space. 

Our input to Eq. ()3.3|) at Y = Yq w 4.6 was taken from the fixed b solution [22], from which we 
evolved over a range of about 10 units of rapidity. The selection of the rapidity difference between 
two successive steps of the evolution was based on the Euler two-step procedure, in which, for 
each i, a first solution is obtained along the path Y^ — > Yi + h, and a second solution is obtained 
along the path Yj, — > Yi + h/2 — > Y± + h. The difference per unit step between the two solutions is 
proportional to h{Y)d 2 N / BY 2 . For each step of the evolution, h(yi) was iteratively reduced until 
the maximal difference between the first and the second solutions over the ith matrix is small (we 
set our accuracy condition to 10~ 2 , and satisfied ourselves that the solution is stable within a few 
percent under variations of this choice). 

The existence of a first derivative of N on both side of Eq. (J3.1)) . slows the convergence rate of the 
equation. Using parallel programming techniques (16 CPUs) and 50 x 50 matrices, an evolution 
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Figure 3: The Y -dependence of the solution to Eq. \3. 1)) for r /R= \ and b = 0, 10 in units of GeV 1 



over 10 units of rapidity was completed within about 50 — 500 hours, depending on the value of 
a 5 (faster convergence for smaller couplings). Due to the numerical complexity of Eq. (|3.1|) . the 
numerical solution presented below is for fixed as- We recognize, however, that QCD evolution 
should be performed with running QCD coupling. 

Fig. El shows the difference between the solution to the pure and modified BK equations. Shown 
in Fig. El are the two solutions, as a function of Y, at r/R — ^ an d two values of the impact 
parameter: 6 = (left) and b = 10 GeV -1 (right). The additional term in Eq. ()3.1|) moderates 
the F-dependence of the solution considerably. The effect is particularly pronounced for large 
b = 10 GeV" 1 . 

Figs. |U and El show the saturation scale, Q sa t(Y,b), as a function of b and Y, respectively. The 
definition of Q sa t is a matter of taste, and we choose to define it as the value of Q = 2/r at which 
21og(l — N(Y,r;b)) = —1. We note that both the value and the slope of Q sa t are significantly 
reduced with the inclusion of the DGLAP term. We have also analyzed the effect of the target 
size on the saturation scale, see Fig. El where we show Q sa t(Y, b) for R 2 = 3, 5 GeV~ 2 . At a 
certain value of the rapidity, which depends on the particular definition of the saturation scale 
and the target size, Q sat exhibits an abrupt transition from a flat, or even decreasing with energy, 
to a steep increasing energy dependence. Such transition was also observed in [3T], in which the 
kernel of the BK equation was modified with exp (— A|r|), where A is a free parameter. Broadly 
speaking, apart from the DGLAP terms and the angle integration, see Eq. (|3.2jl . this approach 
is similar to a more simplified approach we adopt here with A ~ 1/R. 

Perhaps, the most reliable test of the validity of the non-linear evolution approach would be a 
comparison to the F 2 data. However, the energy dependence of the solutions [2H1 EUj to the 
pure BK equation are far too steep to allow such a comparison. The calculation of F 2 involves 
integrating the dipole amplitude over the impact parameter. The increase of F 2 with rapidity can 
therefore be assessed by calculating the F-dependence of (b 2 ), defined as: 



Jd 2 bb 2 N(Y,r;b) 
J d 2 bN(Y, r; b) ' 



(3.4) 



10 



1000 




BK 

BK-modified 



Figure 4: The b- dependence of the saturation scale, b in units of GeV 1 . 
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b=0 




Figure 6: The effect of the target size, R, (in units of GeV 1 ) on the solution of the modified BK 
equation. 
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Fig. shows (b 2 ) as a function of Y for the solutions to the modified and pure BK equations. For 
large Y — 14 the overall effect of the DGLAP term reaches about 60%. 

As stated, QCD evolution should be performed, in principle, for running as- Such a numerical 
procedure is beyond the scope of the present paper. Nevertheless, one can assess the influence 
of the QCD coupling on the evolution process, by performing the evolution with different (fixed) 
values of as- Fig. |S] shows the evolution process at b = 0, for as = 0.2 and as = 0.4. Note that 
both Y- and r-dependences are considerably steeper for large QCD coupling. It is much harder 
to achieve numerical stability for large as- This is due to the DGLAP term in Eq. ([3.1)1 . which is 
proportional to asdN/dY . For small as, the DGLAP term is small compared to the l.h.s. of the 
equation, and the error due to the discretization ([3.3)1 is small. For larger values of as, this error 
becomes greater, resulting in a very small step size h{Yj). The practical consequence of this is an 
increase of about an order of magnitude in computing time for as = 0.4, compared to as = 0.2. 



4 Conclusions 

Our numerical solution to Eq. (J2.14j) . presented in Figs. 0HH1 shows that the influence of the pre- 
asymptotic corrections, related to the full anomalous dimension of the DGLAP equation, is rather 
large. These corrections moderate the energy behaviour of the amplitude, as well as the value 
of the saturation scale. The BK equation without any modification does not provide a reliable 
prediction for the LHC energies, while there is a reasonable chance that the modified equation 
will be more successful. 

We note, however, that we have not addressed two issues: (i) the solution of the full equation 
Eq. ([2.16)1 . and (ii) the inclusion of a running QCD coupling in the numerical procedure. In 
particular, our experience shows that it will be necessary to employ new numerical methods to 
solve Eq. ([2.16)1 . without which, one cannot find a reliable solution. These issues will be the 
subject of a future publication. 
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Figure 8: The evolution process at b = according to the on the modified BK equation, for two fixed 
values of the QCD coupling. 



The running QCD coupling should be used to provide reliable estimates of physical observables 
in the LHC range of energies (see, for example, Refs. |2"2"| EH] and especially Ref. jH2| where the 
running as case is discussed in detail). Fig. |H1 indicates that running as will lead to further 
suppression of the steep increasing energy behaviour (see Ref.jH2] for numerical estimates of this 
suppression) . 

We firmly believe that the modified BK equation, proposed in this paper, can be a basis for a 
new global fit of the experimental data, that will provide reliable predictions for the deep inelastic 
parton densities, in the LHC energy range. 
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